Загрузите датасет very_low_birthweight.RDS
Это данные о 671 младенце с очень низкой массой тела (<1600 грамм),
собранные в Duke University Medical Center доктором Майклом О’Ши c 1981
по 1987 г.
Переменными исхода являются колонки ‘dead’, а также время от рождения до
смерти или выписки (выводятся из ‘birth’ и ‘exit’. 7 пациентов были
выписаны до рождения).
Сделайте копию датасета, в которой удалите колонки с количеством пропусков больше 100, а затем удалите все строки с пропусками.
data <- read_rds("./very_low_birthweight.rds")
data_filtered_100 <- data[, colSums(is.na(data)) <= 100]
#data_no_na <- data_filtered_100[complete.cases(data_filtered_100),]
data_no_na <- data_filtered_100 %>% filter(rowSums(is.na(data_filtered_100)) == 0)
Постройте графики плотности распределения для числовых переменных. Удалите выбросы, если таковые имеются. Преобразуйте категориальные переменные в факторы. Для любых двух числовых переменных раскрасьте график по переменной ‘inout’.
2.1 - Выделяем только числовые переменные и строим графики плотности распределения.
numeric_variables_only <- data_no_na[, sapply(data_no_na, is.numeric)]
for (var in names(numeric_variables_only)) {
plot <- ggplot(numeric_variables_only, aes_string(x = var)) +
geom_density(fill = "blue", alpha = 0.5) +
labs(title = paste("Плотность распределения:", var), x = var) +
theme_minimal()
print(plot)
}
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
2.2 - удаляем выбросы в количественных переменных
Используя метод IQR, я получила очень строгую фильтрацию (из 531 осталось 232 значения, включая всех dead, что важно далее. Я пробовала использовать менее строгие параметры - 2 или 3 вместо 1.5 при умножении), однако это практически не меняет ситуацию. Тут возможно не исключать строку с outlier полностью, а например заменять на NA, но я решила попробовать другой метод.
data_no_outliers <- data_no_na
for (var in names(numeric_variables_only)) {
q1 <- quantile(data_no_outliers[[var]], 0.25, na.rm = TRUE)
q3 <- quantile(data_no_outliers[[var]], 0.75, na.rm = TRUE)
iqr <- q3 - q1
upper <- q3 + 1.5 * iqr
lower <- q1 - 1.5 * iqr
data_no_outliers <- data_no_outliers %>%
filter(data_no_outliers[[var]] >= lower & data_no_outliers[[var]] <= upper)
}
В качестве альтернативы я использовала Z-score и считала outliers только те значения, которые более чем на 3 стандартных отклонения отклоняются от среднего. Тогда из 531 значения дальше проходят 513, что кажется более разумным результатом.
data_no_outliers <- data_no_na
for (var in names(numeric_variables_only)) {
numeric_data <- data_no_outliers[[var]]
z_scores <- scale(numeric_data)
data_no_outliers <- data_no_outliers %>%
filter(abs(z_scores) <= 3)
}
## Warning: Using one column matrices in `filter()` was deprecated in dplyr 1.1.0.
## ℹ Please use one dimensional logical vectors instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
2.3 - переводим категориальные переменные в факторы
categorical_variables <- sapply(data_no_outliers, is.character)
data_no_outliers[, categorical_variables] <- lapply(data_no_outliers[, categorical_variables], factor)
data_no_outliers %>% glimpse() %>% head()
## Rows: 513
## Columns: 19
## $ birth <dbl> 81.514, 81.558, 81.593, 81.610, 81.624, 81.626, 81.684, 81.68…
## $ exit <dbl> 81.539, 81.667, 81.599, 81.697, 81.700, 81.730, 81.854, 81.87…
## $ hospstay <int> 9, 40, 2, 32, 28, 38, 62, 69, 1, 93, 44, 66, 65, 44, 70, 85, …
## $ lowph <dbl> 7.250000, 7.250000, 6.969997, 7.320000, 7.160000, 7.039997, 7…
## $ pltct <int> 244, 182, 54, 282, 153, 229, 182, 361, 378, 255, 186, 260, 18…
## $ race <fct> white, black, black, black, black, white, black, white, white…
## $ bwt <int> 1370, 1480, 925, 1255, 1350, 1310, 1110, 1180, 970, 770, 1490…
## $ gest <dbl> 32.0, 32.0, 28.0, 29.5, 34.0, 32.0, 28.0, 28.0, 28.0, 26.0, 3…
## $ inout <fct> born at Duke, born at Duke, born at Duke, born at Duke, born …
## $ twn <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0…
## $ delivery <fct> abdominal, vaginal, abdominal, vaginal, abdominal, vaginal, v…
## $ apg1 <int> 7, 8, 5, 9, 4, 6, 6, 6, 2, 4, 8, 1, 8, 5, 9, 9, 9, 6, 1, 6, 1…
## $ vent <int> 0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1…
## $ pneumo <int> 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0…
## $ pda <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ cld <int> 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1…
## $ year <dbl> 81.51471, 81.55847, 81.59406, 81.61053, 81.62421, 81.62695, 8…
## $ sex <fct> female, male, female, female, female, male, male, male, femal…
## $ dead <int> 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0…
## birth exit hospstay lowph pltct race bwt gest inout twn
## 1 81.514 81.539 9 7.250000 244 white 1370 32.0 born at Duke 0
## 2 81.558 81.667 40 7.250000 182 black 1480 32.0 born at Duke 0
## 3 81.593 81.599 2 6.969997 54 black 925 28.0 born at Duke 0
## 4 81.610 81.697 32 7.320000 282 black 1255 29.5 born at Duke 0
## 5 81.624 81.700 28 7.160000 153 black 1350 34.0 born at Duke 0
## 6 81.626 81.730 38 7.039997 229 white 1310 32.0 born at Duke 0
## delivery apg1 vent pneumo pda cld year sex dead
## 1 abdominal 7 0 0 0 0 81.51471 female 0
## 2 vaginal 8 0 0 0 0 81.55847 male 0
## 3 abdominal 5 1 1 0 0 81.59406 female 1
## 4 vaginal 9 0 0 0 0 81.61053 female 0
## 5 abdominal 4 0 0 0 0 81.62421 female 0
## 6 vaginal 6 1 0 0 0 81.62695 male 0
2.4 - Визуализация двух числовых переменных с раскраской по inout
ggplot(data_no_outliers, aes_string(x = "birth", y = "hospstay", color = "inout")) +
geom_point(alpha = 0.5) +
labs(title = "Визуализация переменных birth, hospstay с раскраской по inout") +
theme_minimal()
В целом видим отсутствие различий в зависимости от inout.
Проведите тест на сравнение значений колонки ‘lowph’ между группами в переменной inout. Вид статистического теста определите самостоятельно. Визуализируйте результат через библиотеку ‘rstatix’. Как бы вы интерпретировали результат, если бы знали, что более низкое значение lowph ассоциировано с более низкой выживаемостью?
Разбив значения lowph согласно переменной inout, видим, что распределения не нормальные, поэтому будем использовать непараметрический тест, в данном случае тест Манна-Уитни (Wilcoxon).
ggplot(data_no_outliers, aes(x = lowph)) +
geom_histogram(fill = "blue", alpha = 0.5) +
facet_wrap(~ inout) +
theme_minimal()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
wilcox_test_out <- data_no_outliers %>%
wilcox_test(lowph ~ inout) %>%
add_significance()
print(wilcox_test_out)
## # A tibble: 1 × 8
## .y. group1 group2 n1 n2 statistic p p.signif
## <chr> <chr> <chr> <int> <int> <dbl> <dbl> <chr>
## 1 lowph born at Duke transported 433 80 23662 0.000000191 ****
wilcox_test_out <- wilcox_test_out %>%
mutate(y.position = max(data_no_outliers$lowph, na.rm = TRUE) * 1.1)
ggboxplot(data_no_outliers, x = "inout", y = "lowph", fill = "inout") +
stat_pvalue_manual(wilcox_test_out, label = "p", tip.length = 0.01) +
theme_minimal() +
labs(
title = "Сравнение групп lowph в зависимости от переменной inout"
)
В
целом, в группе transported значения ниже, чем в born at Duke. Если
более низкое значение lowph ассоциировано с более низкой выживаемостью,
то возможно, транспортируют более тяжелых больных, либо же сам процесс
транспортировки негативно влияет на параметры пациента.
Сделайте новый датафрейм, в котором оставьте только континуальные или ранговые данные, кроме ‘birth’, ‘year’ и ‘exit’. Сделайте корреляционный анализ этих данных. Постройте два любых типа графиков для визуализации корреляций.
Оставляем континуальные данные и ранговые (apg1) данные, строим корреляционную матрицу.
continuous <- data_no_outliers %>%
select(hospstay, lowph, pltct, bwt, gest, apg1)
cor_matrix <- cor(continuous, method = "pearson", use = "complete.obs")
ggcorrplot(cor_matrix, hc.order = TRUE, type = "lower", lab = TRUE, lab_size = 3,
title = "Корреляции между различными переменными")
Визуализируем корреляции 2 способами - диаграмма рассеяния с линией
регрессии и heatmap.
ggplot(continuous, aes(x = gest, y = bwt)) +
geom_point(color = "blue", alpha = 0.6) +
geom_smooth(method = "lm", color = "red", se = FALSE) +
theme_minimal() +
labs(
title = "Диаграмма разброса с линией регрессии для переменных с максимальной корреляцией в датасете"
)
## `geom_smooth()` using formula = 'y ~ x'
heatmap(scale(continuous))
Постройте иерархическую кластеризацию на этом датафрейме, используя euclidean distances.
continuous_scaled <- scale(continuous)
hierarchical <- hclust(dist(continuous_scaled, method = "euclidean"), method = "ward.D2")
plot(hierarchical, xlab = "", sub = "", labels = FALSE)
Сделайте одновременный график heatmap и иерархической кластеризации. Интерпретируйте результат.
Также используем евклидово расстояние.
pheatmap(cor_matrix, clustering_distance_rows = "euclidean",
clustering_distance_cols = "euclidean", clustering_method = "ward.D2")
Продолжительность пребывания в госпитале обратно корррелирует с остальными показателями. В целом, это можно интерпретировать так, что чем “лучше” показатели, тем меньше требуется находиться в госпитале. Переменные bwt и gest наиболее скоррелированы. Они отражают вес при рождении и срок беременности, так что это вполне естественно.
Проведите PCA анализ на этих данных. Проинтерпретируйте результат. Нужно ли применять шкалирование для этих данных перед проведением PCA?
pca_out <-prcomp(continuous_scaled, center = TRUE, scale. = TRUE)
summary(pca_out)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6
## Standard deviation 1.6049 1.0405 0.8716 0.8610 0.7590 0.51449
## Proportion of Variance 0.4293 0.1804 0.1266 0.1235 0.0960 0.04412
## Cumulative Proportion 0.4293 0.6097 0.7363 0.8599 0.9559 1.00000
summary(pca_out)$importance[2, ]
## PC1 PC2 PC3 PC4 PC5 PC6
## 0.42928 0.18045 0.12661 0.12355 0.09600 0.04412
В этом случае применяем шкалирование, так как данные очень разные, имеют разную размерность. Шкалирование поможет нивелировать влияние параметров с большим абсолютным значением.
Смотрим на процент вариации, который объясняют компоненты. Первая объясняет 43%, первые три вместе - около 73%. Первые две являются наиболее важными.
Постройте biplot график для PCA. Раскрасьте его по значению колонки ‘dead’.
biplot_graph <- cbind(continuous, pca_out$x)
biplot_graph$dead <- data_no_outliers$dead
ggplot(biplot_graph, aes(x = PC1, y = PC2, color = as.factor(dead))) +
geom_point(alpha = 0.5) +
labs(title = "Biplot colored by survival",
x = "PC1", y = "PC2") +
scale_color_manual(values = c("blue", "red")) +
theme_minimal()
Не
видно различий между выжившими и нет.
Переведите последний график в ‘plotly’. При наведении на точку нужно, чтобы отображалось id пациента.
biplot_graph$id <- rownames(biplot_graph)
ggplot_biplot <- ggplot(biplot_graph, aes(x = PC1, y = PC2, color = as.factor(dead), text = id)) +
geom_point(alpha = 0.5) +
labs(title = "Biplot colored by survival",
x = "PC1", y = "PC2") +
scale_color_manual(values = c("blue", "red")) +
theme_minimal()
plotly_biplot <- ggplotly(ggplot_biplot, tooltip = "text")
plotly_biplot
Дайте содержательную интерпретацию PCA анализу. Почему использовать колонку ‘dead’ для выводов об ассоциации с выживаемостью некорректно?
Мы видим, что данные не кластеризуются, PCA не способен в данном случае различать категории выживших и нет. Некорректность выводов по колонке dead может быть связана с тем, что это бинарный показатель и он недостаточно информативен для выводов о выживаемости. Кроме того, мы только знаем, что пациент был жив в какой-то момент времени, но не знаем, что было в дальнейшем.
Приведите ваши данные к размерности в две колонки через UMAP. Сравните результаты отображения точек между алгоритмами PCA и UMAP.
umap_result <- umap(continuous_scaled, n_neighbors = 10, min_dist = 0.1, metric = "euclidean")
biplot_graph_umap <- data.frame(UMAP1 = umap_result$layout[, 1], UMAP2 = umap_result$layout[, 2])
biplot_graph_umap$dead <- data_no_outliers$dead
biplot_graph_umap$id <- rownames(biplot_graph_umap)
ggplot_biplot_umap <- ggplot(biplot_graph_umap, aes(x = UMAP1, y = UMAP2, color = as.factor(dead), text = id)) +
geom_point(alpha = 0.5) +
labs(title = "UMAP Biplot colored by survival",
x = "UMAP1", y = "UMAP2") +
scale_color_manual(values = c("blue", "red")) +
theme_minimal()
plotly_biplot_umap <- ggplotly(ggplot_biplot_umap, tooltip = "text")
plotly_biplot_umap
Как и в случае с PCA, не видим особенных различий между категориями. График выглядит слегка по-другому, имеет несколько другую форму, однако мы не видим закономерностей ни в одном, ни в другом случае.
Давайте самостоятельно увидим, что снижение размерности – это группа методов, славящаяся своей неустойчивостью. Измените основные параметры UMAP (n_neighbors и min_dist) и проанализируйте, как это влияет на результаты.
umap_result <- umap(continuous_scaled, n_neighbors = 50, min_dist = 0.1, metric = "euclidean")
biplot_graph_umap <- data.frame(UMAP1 = umap_result$layout[, 1], UMAP2 = umap_result$layout[, 2])
biplot_graph_umap$dead <- data_no_outliers$dead
biplot_graph_umap$id <- rownames(biplot_graph_umap)
ggplot_biplot_umap <- ggplot(biplot_graph_umap, aes(x = UMAP1, y = UMAP2, color = as.factor(dead), text = id)) +
geom_point(alpha = 0.5) +
labs(title = "UMAP Biplot colored by survival",
x = "UMAP1", y = "UMAP2") +
scale_color_manual(values = c("blue", "red")) +
theme_minimal()
plotly_biplot_umap <- ggplotly(ggplot_biplot_umap, tooltip = "text")
plotly_biplot_umap
При изменении параметров мы видим, что график меняется. Например, теперь представители категории dead теперь скорее сгруппированы в другой его части. Также будто бы паттерн более крупный.
umap_result <- umap(continuous_scaled, n_neighbors = 5, min_dist = 0.1, metric = "euclidean")
biplot_graph_umap <- data.frame(UMAP1 = umap_result$layout[, 1], UMAP2 = umap_result$layout[, 2])
biplot_graph_umap$dead <- data_no_outliers$dead
biplot_graph_umap$id <- rownames(biplot_graph_umap)
ggplot_biplot_umap <- ggplot(biplot_graph_umap, aes(x = UMAP1, y = UMAP2, color = as.factor(dead), text = id)) +
geom_point(alpha = 0.5) +
labs(title = "UMAP Biplot colored by survival",
x = "UMAP1", y = "UMAP2") +
scale_color_manual(values = c("blue", "red")) +
theme_minimal()
plotly_biplot_umap <- ggplotly(ggplot_biplot_umap, tooltip = "text")
plotly_biplot_umap
При уменьшении n_neighbors точки становятся ближе друг к другу, возможно получится различать меньшие различия.
umap_result <- umap(continuous_scaled, n_neighbors = 10, min_dist = 0.5, metric = "euclidean")
biplot_graph_umap <- data.frame(UMAP1 = umap_result$layout[, 1], UMAP2 = umap_result$layout[, 2])
biplot_graph_umap$dead <- data_no_outliers$dead
biplot_graph_umap$id <- rownames(biplot_graph_umap)
ggplot_biplot_umap <- ggplot(biplot_graph_umap, aes(x = UMAP1, y = UMAP2, color = as.factor(dead), text = id)) +
geom_point(alpha = 0.5) +
labs(title = "UMAP Biplot colored by survival",
x = "UMAP1", y = "UMAP2") +
scale_color_manual(values = c("blue", "red")) +
theme_minimal()
plotly_biplot_umap <- ggplotly(ggplot_biplot_umap, tooltip = "text")
plotly_biplot_umap
При увеличении min_dist как будто бы теряется паттерн, точки более рассеянны.